ohibc logo
OHI British Columbia | OHI Science | Citation policy

1 Summary: OHIBC Lasting Special Places subgoal (Sense of Place)

Currently, the Lasting Special Places goal model is identical to the OHI Global model: a region’s status is based upon percent of protected area within 1 km inland buffer and percent of protected area within 3 nautical mile offshore buffer, compared to a reference point of 30% protected area.

\[X_{LSP} = \frac{\frac{pA_{CMPA}}{pA_{refCMPA}} + \frac{pA_{CP}}{pA_{refCP}}}{2}\]

pA = percent of area within the inland or offshore buffer; CMPA = coastal marine protected area (3nm offshore); CP = coastline protected (1km inland); and refCMPA = refCP = 30% reference point for both measures.

Future changes may incorporate other data sets and MaPP planning zones.

2 Data Sources

WDPA data base

3 Methods

3.1 Read in BC WDPA-MPA shapefile

If the WDPA-MPA shapefile has already been rasterized for the global assessment, then we may be able to use that directly. It is Mollweide projection which is not ideal for BC, as we would prefer to work in BC Albers. Some potential solutions:

  • crop the Mollweide to a smaller region enclosing BC EEZ, then reproject it to BC Albers
  • create a region raster for BC in Mollweide (note: already have these in BC Albers)
  • modifying the arcpy script for global LSP, create a BC-specific polygon set, then rasterize to BC Albers.
  • in R, create a BC-specific polygon set, then rasterize to BC Albers.

We’re going with the fourth option.

NOTE: If BC WDPA file does not yet exist, get_wdpa_poly() creates it from the original WDPA-MPA file. This takes a long time, due to reading in the full WDPA-MPA geodatabase into a SpatialPolygonsDataFrame.

poly_wdpa_bc <- get_wdpa_poly(p4s_bcalb, reload = FALSE)  ### defaults to BC Albers

3.2 Rasterize the BC WDPA-MPA shapefile

rast_eez <- raster(file.path(dir_spatial, 'ohibc_rgn_raster_500m.tif'))

wdpa_bc_shp_file  <- file.path(dir_goal_anx, 'int', 'wdpa_bc_bcalb.shp')
wdpa_bc_rast_file <- file.path(dir_rast, 'rast_wdpa_bc.tif')

lsp_rasterize(wdpa_bc_shp_file, 
              wdpa_bc_rast_file, 
              rast_eez, 'STATUS_YR')
## class       : RasterBrick 
## dimensions  : 3126, 3424, 10703424, 1  (nrow, ncol, ncell, nlayers)
## resolution  : 500, 500  (x, y)
## extent      : 159000, 1871000, 173000, 1736000  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## data source : /home/ohara/github/ohibc/prep/lsp/v2016/raster/rast_wdpa_bc_tmp.tif 
## names       : rast_wdpa_bc_tmp
## NULL
rast_wdpa_bc <- raster::raster(wdpa_bc_rast_file)

3.3 Set up coastal buffer rasters

Buffer shapefiles are located in github/ohibc/prep/spatial. LSP uses 1 km inland and 3nm offshore buffers, while resilience requires analysis over the entire EEZ.

Analysis will be done using raster::crosstab() comparing the WDPA raster to various region rasters. Using a 500 m raster is the coarsest that should be used on a 1 km feature; a base raster is available at ~/github/ohibc/prep/spatial/ohibc_rgn_raster_500m.tif.

  • If rasters are not already available for 1 km inland, 3 nm offshore, and EEZ:
    • Read in buffer shapefiles to SpatialPolygonsDataFrames
    • rasterize to same extents/resolution as 500m base raster.
library(tmap)

rast_map <- tm_shape(rast_3nm, is.master = TRUE) +
  tm_raster(alpha = 1, palette = 'Blues') + 
  tm_shape(rast_1km) +
  tm_raster(alpha = 1, palette = 'Greens') + 
  tm_shape(rast_wdpa_bc) +
  tm_raster(alpha = .5, palette = 'Reds')

print(rast_map)

4 Calculate goal model

Once the WDPA raster is cross-tabulated against the OHI region rasters (both 3 nm offshore and 1 km inland) we have the number of protected cells, identified by year of protection, within each region. NA values are unprotected cells.

4.0.1 Summary of zonal stats dataframes (3 nm offshore):

##       year          rgn_id         n_cells       
##  Min.   :1886   Min.   :1.000   Min.   :      0  
##  1st Qu.:1947   1st Qu.:2.000   1st Qu.:      0  
##  Median :1970   Median :4.000   Median :      0  
##  Mean   :1968   Mean   :4.143   Mean   :  15204  
##  3rd Qu.:1992   3rd Qu.:6.000   3rd Qu.:      0  
##  Max.   :2013   Max.   :8.000   Max.   :9942436  
##  NA's   :8      NA's   :88

4.0.2 Summary of zonal stats dataframes (1 km inland):

##       year          rgn_id         n_cells        
##  Min.   :1886   Min.   :1.000   Min.   :       0  
##  1st Qu.:1947   1st Qu.:2.000   1st Qu.:       0  
##  Median :1970   Median :4.000   Median :       0  
##  Mean   :1968   Mean   :4.143   Mean   :   15204  
##  3rd Qu.:1992   3rd Qu.:6.000   3rd Qu.:       3  
##  Max.   :2013   Max.   :8.000   Max.   :10012355  
##  NA's   :8      NA's   :88

4.1 Calculate protected area and total area by region

Grouping by rgn_id, the total number of cells per region is determined by summing cell counts across ALL years, including cells with year == NA (unprotected cells). We can then determine the protected area for each year by looking at the cumulative sum of cells up to any given year.

Since the cells are 500 m on a side, we can easily calculate area by multiplying cell count * 0.25 km2 per cell.

Finally we can calculate the status of a region for any given year by finding the ratio of protected:total and normalizing by the goal’s target of 30% protected area.

4.1.1 Protected areas and status (3 nm offshore, 2010+ only):

4.1.2 Protected areas and status (1 km inland, 2010+ only):

4.2 Combine scores for inland and offshore, and writing output layers

The status is based on a simple arithmetic average of the inland and offshore status values.

We can save outputs for the following layers:

a_prot_inland_file <- file.path(dir_goal, ‘output’, ‘lsp_protected_inland1km.csv’) a_prot_offshore_file <- file.path(dir_goal, ‘output’, ‘lsp_protected_offshore3nm.csv’) a_tot_inland_file <- file.path(dir_goal, ‘output’, ‘lsp_a_total_inland1km.csv’) a_tot_offshore_file <- file.path(dir_goal, ‘output’, ‘lsp_a_total_offshore3nm.csv’)

  • ~/github/ohibc/prep/lsp/v2016/output/lsp_protected_inland1km.csv: inland protected area (km2) for each region (since 1980)
  • ~/github/ohibc/prep/lsp/v2016/output/lsp_protected_offshore3nm.csv: offshore protected area (km2) for each region (since 1980)
  • ~/github/ohibc/prep/lsp/v2016/output/lsp_a_total_inland1km.csv: inland 1 km total area (km2) for each region
  • ~/github/ohibc/prep/lsp/v2016/output/lsp_a_total_offshore3nm.csv: offshore 3 nm total area (km2) for each region

From these layers, we can also estimate the status and trend. “Official” values will be determined in the toolbox? Trend is based on linear model going back ten years from each status year to smooth trend values, since addition of new MPAs is rather sporadic.

Year-by-year status and trend estimates will be saved:

  • ~/github/ohibc/prep/lsp/v2016/output/lsp_status.csv: estimate of status by region since 1980
  • ~/github/ohibc/prep/lsp/v2016/output/lsp_trend.csv: estimate of trend by region since 1990

4.2.1 Status and trend estimates:

4.2.2 Plot map of status by region

Examining OHIBC Lasting Special Places scores for 1995, 2005, and 2015:

4.3 Plot scores time series

To examine results, we plot the estimated status and trend over time.

5 Provenance